An Implicit Time-stepping Method for Quasi-rigid Multibody Systems with Intermittent Contact

We recently developed a time-stepping method for simulating rigid multi-body systems with intermittent contact that is implicit in the geometric information [1]. In this paper, we extend this formulation to quasi-rigid or locally compliant objects, i.e., objects with a rigid core surrounded by a compliant layer, similar to Song et al. [2]. The difference in our compliance model from existing quasi-rigid models is that, based on physical motivations, we assume the compliant layer has a maximum possible normal deflection beyond which it acts as a rigid body. Therefore, we use an extension of the Kelvin-Voigt (i.e. linear springdamper) model for obtaining the normal contact forces by incorporating the thickness of the compliant layer explicitly in the contact model. We use the Kelvin-Voigt model for the tangential forces and assume that the contact forces and moment satisfy an ellipsoidal friction law. We model each object as an intersection of convex inequalities and write the contact constraint as a complementarity constraint between the contact force and a distance function dependent on the closest points and the local deformation of the body. The closest points satisfy a system of nonlinear algebraic equations and the resultant continuous model is a Differential Complementarity Problem (DCP). This enables us to formulate a geometrically implicit time-stepping scheme for solving the DCP which is more accurate than a geometrically explicit scheme. The discrete problem to be solved at each time-step is a mixed nonlinear complementarity problem. INTRODUCTION To automatically plan and execute tasks involving intermittent contact, one must be able to accurately predict the object motions in such systems. Applications include haptic interactions, collaborative human-robot manipulation, such as rearranging the furniture in a house, as well as industrial automation, such as simulation of parts feeders. These applications are characterized by intermittency of contact, presence of stick-slip frictional behavior and deformation at the contact surfaces. The deformation at the contact is usually very small and therefore the objects can be modeled as quasi-rigid or locally compliant [2–5], i.e., each body consists of a rigid core surrounded by a thin compliant shell. Such objects may have a maximum possible deflection and the contact will behave as a rigid contact once the maximum deflection is reached. This motivates us to model the objects as locally compliant objects with a limit on the allowable deflection at the contact. The dynamics of multi-rigid-body systems with unilateral contacts can be modeled as differential algebraic equations (DAE) [6] if the contact interactions (sliding, rolling, or separating) at each contact are known. However, in general, the contact interactions are not known a priori, but rather are discovered as part of the solution process. To handle the many possibilities in a rigorous theoretical and computational framework, the model is formulated as a differential complementarity problem [7, 8]. The primary sources of stability and accuracy problems in dynamic simulation are polyhedral approximations of smooth bodies, the decoupling of collision detection from the solution of the 1 Copyright c © 2007 by ASME dynamic time-stepping subproblem, and approximations to the quadratic Coulomb friction model. Irrespective of whether the model of the object is rigid or locally compliant, all state-of-theart time steppers [4,9,10] are explicit with respect to the geometric information, i.e., they use the geometric information obtained from a collision detection algorithm at the start of the current time step to compute the state at the end of the time step without modifying this information. The method of Tzitzouris [11] is the only geometrically implicit method developed to date, but it requires a closed form distance function between the two bodies which is usually not available. In our previous work [1] we showed simulation results of a disk rolling on a plane without slip and depicted the loss of energy due to polyhedral approximation and the approximation of the distance function. To overcome this, we presented a geometrically implicit time-stepping scheme for convex objects described by implicit surfaces in which the collision detection and dynamic time stepping problem is solved in the same time step. The main focus of this paper is to develop a geometricallyimplicit time-stepping model for dynamic simulation of convex objects described by implicit surfaces, assuming single point contact between the objects and local compliance at the contacts. However, unlike other locally compliant models, we assume a limit on the maximum amount of allowable deflection which is realistic in many scenarios (e.g., to model flesh and bone for biomechanics and human robot interaction). We extend our formulation for contact constraints presented in [1] to include the deflection at the contact and use a linear viscoelastic Kelvin-Voigt model (i.e., a linear spring-damper model) for modeling the compliance. The contact constraints also take into account the maximum allowable deflection at the contact point. Since we are assuming an upper bound on the deflection, there can be an instantaneous jump in the contact forces when the rigid core is reached. Thus we formulate our time-stepping problem at the velocity-impulse level instead of the force-acceleration level so that the resulting time-steppers are better behaved. Our paper is organized as follows. In Section 1, we survey the relevant literature. In Section 2, we present both the continuous and discrete time dynamics model for multi-rigid-body systems with an ellipsoidal dry friction law. In Section 3, we review the non-penetration condition for the contact constraints presented in [1]. Thereafter, in section 4, we modify these contact constraints to include compliant contacts with limits on the maximum allowable deflection. The discrete time dynamics model along with the contact constraints form a mixed nonlinear complementarity problem at each time-step. In Section 5, we give examples that validate and elucidate our time-stepping scheme. Finally in section 6, we present our conclusions and lay out the future work. RELATED WORK Dynamics of multi-rigid-body systems with unilateral contacts can be modeled as differential algebraic equations (DAE) [6] if the contact interactions (sliding, rolling, or separating) at each contact are known. However, in general, the contact interactions are not known a priori, but rather are discovered as part of the solution process. To handle the many possibilities in a rigorous theoretical and computational framework, the model is formulated as a differential complementarity problem [7, 8]. The differential complementarity problem is solved using a timestepping scheme and the resultant system of equations to be solved at each step is a (linear/nonlinear) complementarity problem. Definition 1 (Nonlinear Complementarity Problem (NCP)). Let f (z) ∈ Rn be a given vector function of z ∈ Rn. The nonlinear complementarity problem is to find z satisfying 0 ≤ z ⊥ f (z) ≥ 0, where the symbol ⊥ connotes orthogonality (i.e., f (z) · z = 0) and vector inequalities hold on a per element basis. When f (z) is linear in z, then the problem is referred to as a linear complementarity problem (LCP). Of particular importance to this work is a generalization of the NCP known as the mixed nonlinear complementarity problem [12]. Definition 2. Let g(u,v) : Rn1 ×Rn2 → Rn1 and f (u,v) : Rn1 × R n2 →Rn2 be given vector functions of u∈Rn1 and v∈Rn2 , with n1 +n2 = n. The mixed nonlinear complementarity problem is to find u and v satisfying g(u,v) = 0, u, free 0 ≤ v ⊥ f (u,v) ≥ 0 Frictional collisions between rigid bodies have a long history in mechanics [13, 14]. Here, we give an overview of the basic approaches and refer the reader to a recent survey article [15] for a more comprehensive review. There are two primary approaches to modeling collisions: coefficient of restitution based approaches and force based methods. In the former, the process of energy transfer and dissipation during collision is modeled by various coefficients relating the velocity (or impulses) before contact to that after contact. However, the extension of these concepts to situations with multiple contacts is not straightforward. The force based approaches use a compliant contact model to compute the contact forces where the contact compliance is modeled as a (linear/nonlinear) spring-damper system. In the simplest model (known as Kelvin-Voigt model or linear spring-damper model), the normal contact force is given by a linear function of the deformation and the rate of deformation (F = kδ + cδ̇) i.e. the flexibility of the body is lumped as a linear spring (with spring constant k) and damper (with damping coefficient c). The limitations of the linear model are documented in [15]. Hertz introduced a nonlinear model of the form 2 Copyright c © 2007 by ASME F = kδn, where n is a constant [16]. This model was extended to a nonlinear spring-damper model by Hunt and Crossley [17] of the form F = kδn + cδpδ̇q, where p,q are constants. The models presented above are believed to be of increasing accuracy but there are more unknown constants dependent on geometry of the objects and material properties that have to be determined experimentally (except for some simple cases). This is a general feature of all proposed contact compliance models. In [18] a continuum model of the deformations at each contact is used. Song and Kumar [2] have used a 3D linear distributed contact model to compute the contact forces. In this paper we use a lumped 3D linear spring-damper to model the contact compliance similar to [19]. However, we note that we could have replaced this with a lumped nonlinear model if required. We use an elliptic dry friction law [20] that is a generalization of Coulomb’s friction law to model the friction at the contact. DYNAMIC MODEL FOR RIGID BODY SYSTEMS In complementarity methods, the instantaneous equations of motion of a multi-rigid-body system consist of five parts: (a) Newton-Euler equations, (b) Kinematic map relating the generalized velocities to the linear and angular velocities, (c) Equality constraints to model joints, (d) Normal contact condition to model intermittent contact, and (e) Friction law. Parts (a) and (b) form a system of ordinary differential equations, (c) is a system of (nonlinear) algebraic equations, (d) is given by a system of complementarity constraints, and (e) can be written as a system of complementarity constraints for Coulomb friction law using the maximum work dissipation principle. In this paper we use a more general elliptic dry friction law [20]. Thus, the dynamic model is a differential complementarity problem (DCP). To solve this system of equations, we set up a time-stepping scheme and solve a complementarity problem at each time step. We present below the instantaneous-time formulation as well as an Euler time-stepping scheme. To simplify the exposition, we ignore the presence of joints or bilateral constraints in the following discussion. However, all of the discussion below holds in the presence of bilateral constraints. To describe the dynamic model mathematically, we first introduce some notation. Let q j be the position and orientation of body j in an inertial frame and ν j be the concatenated vector of linear (v) and angular (ω) velocities. The generalized coordinates, q, and generalized velocity, ν of the whole system are formed by concatenating q j and ν j respectively. Let λin be the normal contact force at the ith contact and λn be the concatenated vector of the normal contact forces. Let λit and λio be the orthogonal components of the friction force on the tangential plane at the ith contact and λt, λo be the respective concatenated vectors. Let λir be the frictional moment about the ith contact normal and λr be the concatenated vector of the frictional moments. Let nb be the number of bodies and nc be the number of contacts. The instantaneous dynamic model can then be written as follows: Newton-Euler Equations of Motion: M(q)ν̇ = Wnλn + Wtλt + Woλo + Wrλr + λapp + λvp (1) where M(q) is the inertia tensor, λapp is the vector of external forces, λvp is the vector of Coriolis and centripetal forces, Wn, Wt, Wo, and Wr are dependent on q and map the normal contact forces, frictional contact forces, and frictional moments to the body reference frame.